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We investigate the impact of nonlinear evolution of the gravitational potentials in the LCDM 
model on the Integrated Sachs-Wolfe (ISW) contribution to the CMB temperature power spectrum, 
and on the cross-power spectrum of the CMB and a set of biased tracers of the mass. We use an 
ensemble of A^-body simulations to directly follow the potentials and compare the results to analytic 
perturbation theory (PT) methods. The predictions from the PT match the results to high precision 
for k < 0.2 /iMpc -1 . We compute the nonlinear corrections to the angular power spectrum and find 
them to be < 10% of linear theory for I < 100. These corrections are swamped by the cosmic 
variance. On scales I > 100 the departures are more significant, however the CMB signal is more 
than a factor 10 3 larger at this scale. Nonlinear ISW effects therefore play no role in shaping the CMB 
power spectrum for I < 1500. We analyze the CMB-density tracer cross-spectrum using simulations 
and renormalized bias PT, and find good agreement. The usual assumption is that nonlinear 
evolution enhances the growth of structure and counteracts the linear ISW on small scales, leading 
to a change in sign of the CMB-LSS cross-spectrum at small scales. However, PT analysis suggests 
that this trend reverses at late times when the logarithmic growth rate / = dlnD/dlna < 0.5 or 
0. m {z) < 0.3. Numerical results confirm these expectations and we find no sign change in ISW-LSS 
cross-power for low redshifts. Corrections due to nonlinearity and scale dependence of the bias 
are found to be < 10% for I < 100, therefore below the signal-to-noise of the current and future 
measurements. Finally, we estimate the cross-correlation coefficient between the CMB and halos 
and show that it can be made to match that for the dark matter and CMB to within 5% for thin 
redshift shells, thus mitigating the need to model bias evolution. 



I. INTRODUCTION 

Measurements of the temperature fluctuations in the 
Cosmic Microwave Background (CMB), provide a unique 
window onto the primordial Universe and a means to 
learn about the physical processes that generated the ini- 
tial conditions. This discriminatory power is exemplified 
by recent results from the WMAP experiment [lj: the 
primordial power spectral index is n s = 0.960 ± 0.013, 
ruling out the Harrison-Zel'Dovich spectrum at 3ct level. 
However, the temperature power spectrum does not pro- 
vide a pristine window, but it must be cleaned for the im- 
print of foreground signals. One cosmological foreground, 
is the change in energy that a CMB photon experiences as 
it propagates through an inhomogeneous Universe with 
time evolving gravitational potentials, There are three 
main effects that may give rise to such secondary fluctu- 
ations: 
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• Linear Integrated Sachs- Wolfe Effect 0, hereafter 
ISW]: unless the growth of density perturbations 
matches the expansion rate, $ will evolve from zero. 
This will lead to a net change in photon tempera- 
tures. In LCDM |<E>| < as the potential decays, 
giving rise to a net positive correlation between 
density and temperature in Fourier space. 

• Rees-Sciama Effect [!, hereafter RS] : nonlinear col- 
lapse of perturbations to filaments and clusters 
leads to 4> ^ even in the absence of linear ISW, 
and CMB photons change energy as they transit 
across nonlinear structures. It is usually assumed 
that nonlinear evolution accelerates the growth of 
structure and counteracts the linear decay of gravi- 
tational potential in LCDM. In this paper we show 
that this is not always justified. 

• Birkinshaw-Gull Effect 0, hereafter BG]: if a mass 
concentration moves transversely to the line of 
sight, it will create a time variation in the potential 
even if the potential itself is not evolving in time, 
and this will have a dipolar pattern. Consequently, 
photons which enter the potential in the wake, will 
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receive a net energy boost on exit, and those which 
enter ahead will loose energy on exit. However, un- 
like the previous two effects, this contributes only 
to the CMB auto-correlation and not to the cross- 
correlation of CMB with a density tracer. 

All three effects combine into the nonlinear 1SW. It is well 
known that the linear ISW effect leads to fluctuations of 
the order AT « 1/iK on the largest scales I < 10 for 
LCDM [see for example H[ and has been used to rule out 
the self-accelerating branch of DGP model • 

The impact of the nonlinear evolution of $ on the 
CMB has been the subject of a number of studies. 
However, most of these works attempt to quantify the 
effect thr oug h the use of simplified analytic models 
[I 0, S, H d d • A number of studies have employed 
numerical simulations to track the evolution of <!>: in a pi- 
oneering study, Tuluie and Laguna fl3j and Tuluie et al. 
[l4| used ray tracing methods to compute the change in 
temperature for individual photon bundles propagating 
through inhomogeneous universes. They found that the 
combined imprint on the CMB power spectrum, due to 
the RS and BG effects, were of the order AT ~ 1/iK 
on angular scales I ~ 200. Owing to the limited size of 
their simulations, they were unable to comment on the 
effects on the lower multipoles. Seljak related 4> to 
density and momentum using the Poisson and continuity 
equation. These predictions were compared to an N- 
body simulation of the then favored SCDM model, and 
good agreement was found between the two as well as to 
those of Tuluie et al. However, these results were 

obtained in the context of Q m = 1 model, where no linear 
ISW exists, and they could not address I < 100 behavior, 
owing to the limited dynamic range of the simulations. 
Puchades et al. [l6| also recently addressed this problem, 
but again attention was focused on the large multipole 
regime. 

In a more recent study, Cai et al. [TtJ used a single 
A-body simulation, the L-BASIC simulation, which has 
N = 488 3 and comoving length of L = 1.34/i _1 Gpc, to 
compute the nonlinear ISW effect. They measured the 
<j> power spectra at each epoch in the simulation and de- 
veloped an empirical fitting formula for the deviations 
from linear theory. Using this model they computed 
the CMB angular power spectrum and found, on scales 
I > 50, that there was significant nonlinear amplifica- 
tion of power, qualitatively confirming the earlier halo 
model predictions of Cooray [ll[. However, these nonlin- 
ear corrections occur on angular scales where the primary 
anisotropy spectrum is more than two orders of magni- 
tude larger, rendering them of negligible importance. Cai 
et al. also found that there was no evidence for devia- 
tions for multipoles I < 50. One of the aims of this paper 
is to place more precise constraints on the expected level 
of contamination on these large scales. 

The temperature fluctuations induced through the 
evolving 4> can also be observed by correlating the CMB 
against density perturbations, as pointed out by Crit- 
tenden and Turok [la ], and the large-scale ISW effect 



provides an important test for Dark Energy and the cur- 
vature of the Universe. This information can be extracted 
through the cross-correlation of the CMB with tracers of 
the Large-Scale Structure (hereafter LSS). This analysis 
has recently been performed by a number of groups using 
the WMAP data and several large-scale structure mea- 
surements (e.g. SDSS, NVSS, 2MASS). This work has 
resulted in up to 4a level detections of the ISW effect 
[3 HI M, El IM [H, HE [H, i3. In the near future 
these detections will be improved upon with PLANCK 
and the new wide field LSS surveys, such as BOSS, DES, 
Pan-STARRS-1 and EUCLID, etc.. However, in a recent 
paper Granett et al. [28| measured the cross-correlation 
between superstructures and super- voids with the CMB. 
On stacking the signal they found a ~ 4.5<t detection, in 
multiple WMAP bands, and the sign of which appeared 
consistent with late time ISW. This appears in stark con- 
trast to expectations from simple signal-to-noise calcula- 
tions within the LCDM model [3 1^ IM SI • A follow 
up 'consistency' test was performed by Granett et al. 
[33 |. the results of which cast some doubt on the the 
signal as arising from ISW, at least within the LCDM 
model. Cai et al. [I?) also investigated the ISW-density 
cross-correlations, focusing on the nonlinearities arising 
from the mass evolution. They found that there was no 
evidence for enhancement of evolution of in agreement 
with the earlier work of Verde and Spergel [33[ . One of 
the questions we shall address in this paper is whether 
selecting biased tracers of LSS relative to the mass dis- 
tribution can influence the detection sensitivity for the 
ISW. 

We pursue a two-pronged attack on all of these prob- 
lems. Our first avenue will be to use a large ensemble 
of A-body simulations to directly follow the evolution 
of <I>. Our second line is analytic, and we use the non- 
linear gravitational perturbation theory (PT) and renor- 
malized bias frameworks to compute all measured quan- 
tities. This will help us to provide physical insight into 
the results along the way. 

The paper breaks down as follows: In fJTT] we summa- 
rize the basic theory of the ISW. In mil] we describe the 
ensemble of simulations that we use, and describe our 
estimator for measuring $ from the simulations. Here 
we also present maps, comparing the time evolution of 
density, and $ in the simulations. In illVI we investigate 
the two-point statistics of and besides the usual linear 
analysis we derive nonlinear expressions within the con- 
text of the gravitational perturbation theory. We evalu- 
ate the theory and compare directly with measurements 
from the simulations. Then in SjV] we compute the im- 
pact on the CMB temperature power spectrum. In WII 
we turn to the cross-correlations with dark matter, fol- 
lowed by the correlations with halos in Will including 
the effects of scale dependent bias. Again, the theory is 
compared directly with measurements from the simula- 
tions. In $ Villi we perform the line-of-sight integrals and 
compute angular cross-power spectra. Finally, in i jlXI we 
summarize our findings and conclude. 
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II. THEORETICAL BACKGROUND 
A. The ISW effect 

On arrival at the observer, the CMB photons, which 
are sourced at the surface of last scattering, z « 1100, 
are imprinted with two sets of fluctuations: the primary 
anisotropies, which are induced by the primordial fluctu- 
ations, perhaps seeded through the inflationary mecha- 
nism; and the secondary anisotropies, which are induced 
as the photons propagate through the clumpy Universe. 
The primary anisotropies have been studied in great de- 
tail for several decades [and for a review of the important 
processes see [H, [Hj]. There are a number of physical 
mechanisms that give rise to the generation of secondary 
anisotropies [for a review see [35[ and one of these is the 
redshifting of the photons as they pass through evolving 
gravitational potentials. 

The temperature fluctuation induced by the gravita- 
tional redshift may be written as Q: 



AT(n) 



dt<f>(n, X ;t) , 



(1) 



where n is a unit direction vector on the sphere, $ is 
the dimensionless metric perturbation in the Newtonian 
gauge, which reduces to the usual gravitational potential 
on small scales, the 'over dot' denotes a partial derivative 
with respect to the coordinate time t from the FLRW 
metric, \ is the comoving radial geodesic distance x = 
J cdt/a(t), and so may equivalently parameterize time, to 
and ti s denote the time at which the photons are received 
and emitted (i.e. last scattering), respectively, c is the 
speed of light and a(t) is the dimensionless scale factor. 

On scales smaller than the horizon, the perturbed Pois- 
son equation enables us to relate potential and matter 
fluctuations 13611 : 



V 2 $(x;i) = 4irGp(t)6(x;t)a 2 (t) 



(2) 



where p(t) is the mean matter density in the Universe and 
the density fluctuation is <5(x;i) = [p(x, t) — p(t)]/p(t). 
Poisson's equation may most easily be solved in Fourier 
space, upon which we have, 



$(k;t) = -4%Gp(t)a 2 (t) 



5(k;t) 



(3) 



However, what we are really interested in is the instan- 
taneous time rate of change of the potential, 



$(k;t) 



AnG 



[p(t)a 3 



d_ 

dt 

H(t) 
a(t) 



*(k;t) 



a(t) 
S(k;t) 



S(k;t) 
a(t) 



(4) 
,(5) 



where [a 3 (t)p(t)] is a time independent quantity in the 
matter dominated epoch. In the above, we also de- 
fined Hit) = a{t)/a{t) and ft m (i) = p(t) / p CI [t(t) , with 



Pcrit(i) = 3H 2 (t) /8ttG. All quantities with a subscript 
are to be evaluated at the present epoch. Estimating 
the change in the photon temperature due to the evolv- 
ing potentials requires knowledge of the evolution of the 
density perturbation and its time rate of change. In the 
linear regime we may solve the equation of motion for 
6 exactly and obtain both of these quantities. However, 
in the nonlinear regime the situation is more complex 
and requires numerical simulations or nonlinear models 
to proceed. In simulations, measuring <5(k, a) is relatively 
straightforward, whereas its time derivative is more com- 
plicated. As was shown by Seljak (l5j one may obtain this 
from the perturbed continuity equation [36| : 



V • [1 + <5(x; t)} v p (x; t) = -a(t)S(x; t) , 



(6) 



where v p (x;i) is the proper peculiar velocity field. On 
defining the pseudo-peculiar momentum field to be, 



p(x; t) = [l + <S(x; f)] v p (x; t) , 



(7) 



then in Fourier space we may solve the continuity equa- 
tion directly to give us: 



S(k;t) = ik-p(k;i)/a(f) 
Hence, our final expression becomes, 



$(k;t) =T(k) 



H{t) £k.p(M) 



a(t) 



a?{t) 



(8) 



(9) 



where to enable us to pass easily from potential to density 
we introduced the quantity 



3 n ( Ho 



(10) 



III. THE ISW FROM iV-BODY SIMULATIONS 

A. The zHQRIZQN simulations 

In this study we use a subset of the Zurich Hori- 
zon, "zHORIZON", simulations. These are a large en- 
semble of pure cold dark matter ./V-body simulations 
(-^Vsim = 30), performed at the University of Zurich on 
the zB0X2 and zB0X3 super-computers. The specific aim 
for these simulations is to provide high precision mea- 
surements of cosmic structures on scales of the order 
~ 100 ft, _1 Mpc and to also provide insight into the rarest 
fluctuations within the LCDM model that we should 
expect to find within the observable universe. In this 
paper we shall only employ the first 8 zHORIZON sim- 
ulations, since these runs have 11 snapshots logarith- 
mically spaced in the expansion factor from z = 1 to 
0, thus giving sufficient time sampling of the simulated 
density field to capture the late time evolution. The 
expansion factors at which snapshots are recorded are: 
a = {1.0, 0.93, 0.87, 0.76, 0.66, 0.62, 0.57, 0.54, 0.5}. 
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Each numerical simulation was performed using the 
publicly available Gadget-2 code [37|, and followed the 
nonlinear evolution under gravity of N — 750 3 equal 
mass particles in a comoving cube of length L = 
1500 /i _1 Mpc. All of the simulations were run within 
the same cosmological model, and the particular choice 
for the parameters was inspired by results from the 
WMAP experiment [J, [H, |3{|. The parameters are: 
{n m0 = 0.25, fb E , = 0.75,ft b ,o = 0.04, a 8 = 0.8, n s = 
1.0, Wo = — 1, h = 0.72}, where these are: the density pa- 
rameters in matter, dark energy and baryons; the power 
spectrum normalization and primordial spectral index; 
equation of state parameter for dark energy p/ p = wq; 
dimensionless Hubble parameter. The transfer function 
for the simulations w as g enerated using the publicly avail- 
able cmbfast code [40j, |4l|, with high sampling of the 
spatial frequencies on large scales. Initial conditions were 
lain down at redshift z = 50 using the serial version of 
the publicly available 2LPT code [43.1431]. 

Dark matter halo catalogs were generated for all snap- 
shots of each simulation using the Friends-of-Friends 
(FoF) algorithm (44[, with the linking- length parame- 
ter set to the standard b — 0.2. For this we used the 
fast parallel B-FoF code, kindly provided by V. Springel. 
The minimum number of particles for which an ob- 
ject was considered to be a bound halo, was set to 
30 particles. This gave a minimum host halo mass of 
- 1.5 x 10 13 M Q //i. 



sity fluctuation is 



N 



NVi 



Nr. 



N 



N 



(13) 



where iV co n = V^/Vw is the total number of grid cells. 

The pseudo-momentum field may be estimated in a 
similar fashion. For convenience we write, 



[1 + S(x)] u(x)a(t) 



(14) 



where u = v p /a is the comoving peculiar velocity field. 
The particle momentum field is then written as 



[(l + <5)u](x) = ^^^(x-x0ui 



N 



N 



(15) 



This may be convolved with the mass assignment scheme 
to obtain the mesh averaged quantity 



[(l + <5)u](x ij . fe ) = i^^u^(x ij)t -x i ) . (16) 

Thus our estimate for the pseudo-momentum field is 
given by 



B. Estimating the ISW effect in simulations 



p(xy fe ) = a(t)^p- ^ uiW{x ijk - x/) 



iV 



N 



(17) 



In order to estimate $, we require estimates of both 
the density field and pseudo-peculiar momentum field in 
Fourier space (c.f. Eq. ((9]). The dark matter density field 
can be written as a sum over Dirac delta functions, 



E 

i=i 



p(x) = 2]mi6 D (x - xi) 



(11) 



where m; is the mass of the Ith particle and we take all 
particles to have equal mass. The density field averaged 
on a cubical lattice can then be obtained through the 
convolution, 

p g (x ijk ) = -!- J d 3 xp(x)W(x i: j k -x) ; 



V, 



^2W{xi jk -xj) 



(12) 



The density Fourier modes were then estimated using 
the publicly available FFTW routines [46], and each re- 
sulting mode was corrected for the convolution with the 
mass-assignment window function. For the CIC algo- 
rithm this corresponds to the following operation: 



5 d (k) =S g (k)/Wcic(k) 



where 



Wcic(k) = n 



i=l,3 



sin [7rfci/2/cN y ] 
[7rfc 4 /2fc; 



NyJ 



(18) 



(19) 



^ i 



and where sub-script d and g denote discrete and grid 
quantities, and where kjq y = irN g /L is the Nyquist fre- 
quency, and N g is the number of grid cells [4 51 ] . 

To obtain the real space 4>(x, t), we solved for $(k, t) 
in Fourier space using Eq. (JOj) , set the unobservable k = 
mode to zero, and inverse transformed back to real space. 



where W represents the dimensionless window function 
of the mass assignment scheme, and where the normaliza- 
tion factor is Vw = J d 3 x'W(x — x'). The filter function 
W that we adopt throughout is the 'cloud-in-cell' charge 
assignment scheme [45| . Hence, our estimate for the den- 



C. Visual representation of the evolution of <& 

Fig.[T]shows how the dark matter particle number, pro- 
jected in a slab of thickness Ax = 100/i _1 Mpc and side 
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FIG. 1: Evolution of 8 in a slab of thickness Ax = 100 h 1 Mpc. The panels, going from left to right and top to bottom, 
represent redshifts: z — {15, 10, 5, 3, 1, 0}. 
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FIG. 2: Evolution of $ in a slab of thickness Ax = 100 h 1 Mpc. The panels, going from left to right and top to bottom, 
represent redshifts: z — {15, 10, 5, 3, 1, 0}. 
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length L — 1500 ft. _1 Mpc, evolves as a function of cosmic 
time from z = 15 to the present day. At early times, one 
can see that the Universe is regular and homogeneous, 
and the imprint of the initial grid configuration is still 
noticeable. At later times, gravitational instability of 
the matter has led to the formation of a pattern of web 
like structures with dense clumps at the vertices of the 
web - the 'Cosmic Web'. The point we wish to stress, is 
that it is difficult for the eye to pick out features that are 
larger than 100 ft, Mpc. 

Fig. [5] shows the evolution of <E>(x, i) as a function of 
cosmic time. At early times [z ~ 15), when there is 
no linear ISW, the maps are dominated by a small-scale 
foam-like structure. At later times, z ~ 10, the foam is 
sharpened and transformed, with butterfly like features 
present at the high density regions, as expected from the 
BG effect, i.e. a flow of mass moving transversely across 
the sky. At later times the dominating structures are on 
extremely large scales (r > 500 /i _1 Mpc), as expected by 
the linear late time ISW effect. 



IV. DENSITY, MOMENTUM AND POTENTIAL 
POWER SPECTRA IN PT 



A. The 3D power spectra 

The perturbed fields of interest may be written as 
Fourier series, 

3 

ipryikj) = ^- y*d 3 a;V' 7 (x)exp[ik J • x] , (21) 

where ip 7 = {<5(x), V • p(x), <fe(x)}, and where is some 
large region of the Universe over which we shall assume 
that the functions obey harmonic boundary conditions. 
Then, from translational invariance and isotropy, the cor- 
relation of different Fourier modes can be written 

P 7172 (fc 4 )C, = V» (k^ 72 (k,)) , (22) 

where -P 7l72 is the power spectrum matrix of all of the 
fields. Using Eq. |9]) we find, for example: 



V§ Pss{k) 



P si (k,a)= T{k) 



H(a) 



Ps b s(k) 



1 



-PsUk) 



(23) 



(24) 



a(i)""" a 2 (t) 
where we have defined u>(k; t) = ik • p(k; t) — 5(k; t)a(t) 



B. Linear theory results 

The two-point statistics may be evaluated easily within 
the linear theory: 5 <C 1 and V ■ v < 1. In this limit 
the Fourier mode of the density and its time derivative 
evolve as: 

*(k;t) = D(t)S(k;t ) ; (25) 
j(k;t) - f{t)H(t)D(t)6(k;t ) , (26) 

where we have the usual definition of the logarithmic 
derivative of the growth factor, 



f(a) =/(O m (a),0 DE (a)) = 



d\ogD(t) 
dloga(t) 



(27) 



Hence we have 
and 



a(t)D(t)\ (\5(k;to)r)V» ; 

a(t)f(a)H(t)D(t)] 2 PM n (k;t Q ) 
a(t)f(a)H(t)] 2 P^(k;t) 



(28) 



Pff(k,t) = o(t)i?(t)D(t)^|*(k;to)| )V M ; 

= {a(t)f(a)H(t)]Pti n (k;t) . (29) 
Inserting these expressions into Eqs (|23[) and gives: 
'H(a) 



Pk 1 £(k,t)=[f(k)Y 



(l-/(o)) 



Ptt n (k;t) (30) 



Pl?(k,t) = f(k) 



H(a) 



(l-/(a)) 



Pjf(M). (31) 



At this point we may note the well known result that, if 
the Universe is in an Einstein-de Sitter (EdS) phase of 
expansion (i.e. O m (a) = 1, and D(a) oc a), then /(a) = 1 
and the bracketed terms in Eqs ([30]) and (|3~lj) vanish, 
so the ISW effect vanishes. However, if the Universe is 
under/overdense in gravitationally active matter, then 
we expect a non-zero signal, which is positive for both 
spectra. In the currently favored LCDM model, Q m o ~ 
0.25, and so 1 — f(a) < 1. However, at early times fi TO — > 
1 and the ISW is shut off. In the next section we explore 
how this picture changes as the fields are evolved into the 
mildly nonlinear regime. 

Before moving on though, we point out that in the lit- 
erature there are a number of commonly used approxima- 
tions for f(a): for example, f(a) w tt 6 [36| : and some- 
what better, f(a) » ft m (a) a6 + [l + ifi m (a)] 
for models with a cosmological constant A; and better 
still the previous formula, but with the power-index of the 
first term 0.6 -> 4/7 (HSi]. However, all these approxi- 
mations deviate at the few percent level when compared 
to the exact result obtained from numerically solving the 
differential equation for linear growth [for further details 
see for example HH, [H(3] ■ We therefore adopt the exact 
numerical solutions for both D(a) and f(a) throughout 
this study. 
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C. Beyond linear theory: Nonlinear PT 

The collapse of cosmic structures can be followed into 
the nonlinear regime using standard perturbation theory 
(PT) methods, applied to an ideal fluid in a uniformly 
expanding spacetime [for an excellent review see[5l|. The 
first application of PT methods to estimate the impact of 
the nonlinear evolution of $ on the CMB, was given by 
Seljak fl5| . That work was conducted within the context 
of the flat SCDM model, and hence only provided an 
estimate for the Rees-Sciama contribution. Furthermore, 
owing to the fact that (1 — /) = at all times in the EdS 
model, it was necessary only to calculate the PT up to 
2nd order in <5, whereas in more general cosmologies, to 
be consistent at first order, one requires the corrections 
up to 3rd order. We shall now calculate the nonlinear 
ISW in the PT framework for the LCDM model. 

To begin, we require from the PT theory the solutions 
for the fluid overdensity, and in Fourier space these may 
be written as, 



*(k,t) = £[D(t)]%,(Mo) 



where the perturbative solutions at each order can be 
written 



<5„(k) 



(27T) 



3n— 3 



qi,-,q* 



(33) 

In the above expression <5i(qj) represents an initial field 
at wavenumber q^, and the nth order perturbed density 
depends on n initial fields. The quantities Fn^ (qi, q„) 
represent the standard PT interaction kernels, sym- 
metrized in all of their arguments. Also we have adopted 
the short-hand notation [<5" D (k)] n = 5 D (k— q^ — ■ ■ — q n ). 
The Dirac delta function ensures that the waves conserve 
momenta through the interaction, i.e. k = qi + • • ■ + q„. 
For example, the second order PT kernel can be written, 



^(qi.qa) 



i 



q-2 



<12 

9i 



7 



Ml,2 



(34) 



where /ii )2 = qi • qY 

In the standard approach of nonlinear PT, the fluid 
equations are solved for the flat EdS background model. 
In this case the spatial and temporal parts of the evolu- 
tion are fully separable and the perturbative solutions at 
each order simply scale as powers of the expansion factor 
a(t) [12]. However, this is not the case for more general 
cosmological models, nevertheless a very good approxi- 
mation to the evolution can be obtained by replacing the 
powers of a(t) — > D(t). Strictly speaking, the PT interac- 
tion kernels also inherit some time dependence, however 
this is an extremely weak function of time and so to a 
very good approximation we may use the kernels from 
the flat EdS case [for deeper discussion of this see [5l| . 

As Seljak [15J showed, to calculate the ISW we simply 
require the PT expansion for S and its time derivative. 



TABLE I: Sign of the NLO correction to the P s ^(k) power 
spectrum. Recall that a positive correction means an increase 
in the decay rate of the potentials. See text for further details. 



Sign 


of correction 


\Pis(k)\ > P 22 (fc) \Pa{k)\ < P 22 (fc) 




> f2 m (a R s) 

< fim(aRs) 


(+) (-) 
(-) (+) 



Using Eq. (|32p . this latter quantity may be written, 

oo 

S(k,t) ^ f(a)H(a)Y,n[D(t)} n 6 n (k : t ) . (35) 

n=l 

These quantities may now be used to compute the Next- 
to-Leading-Order (NLO) corrections to the power spec- 
tra. Using the above expressions plus the standard PT 
techniques we find, for pseudo-momentum: 



pNL pLin _i_ plLoop . 

LOLO LOLO ' LOLO ' 



pNL 

r &LO 



pLin I plLoop , 



(36) 
(37) 



(32) where the one-loop corrections are, 



^°° P = {af(a)H{a)] 2 [APll{k,a)+ZP}l{k,a)]{^) 
P£°° P = laf(a)H(a)]2P£°° p (k,a) . (39) 



For the <f> we find: 



pNL 

ii 

pNL 

si 



pLin 



P 



Lin 
<54> 



plLoop _ 
plLoop 

si ' 



(40) 
(41) 



where the one-loop corrections are, 



plLoop 



(k) 



{[l-2/(a)] 2 P^(fc) 



P 5 lLoo P(fc) 



+ [l-3/(a)][l-/(a)]P 5 1 |(fc)} ; (42) 
T(k)^{l-2f(a)}p£°°v(k,al43) 



In the above expression P 5 2 | and P$$ are the NLO cor- 
rections to the matter power spectrum and we defined, 



P 



lLoop 



,1A 



(k,a) 



pl3 

r ss 



(k,a) 



fc , a) [fo r explicit forms 



for the lLoop expressions, seel52l [53, \54\. 

We may now learn how the NLO corrections entangle 
the pure linear ISW decay of potentials with the nonlin- 
ear RS effects. The easiest way to discern the changes is 
to consider the sign of the corrections in the above equa- 
tions. We notice that there are two ways the sign may 
change: firstly, there is a sign flip with scale, since Pgg 
is negative and dominant on large scales and Pgg is pos- 
itive and dominant on smaller scales; secondly, the time 
dependent prefactors may change sign. 

Considering the time dependent factors, we see that 
the cross-power spectrum, Eq. f|43[) , will only change sign 
when, 1 — 2/(a) = 0, which occurs when f2 « 0.3. We 
shall label this time ors- Table[T]summarizes the changes. 
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The key point to notice is that at early times, there is 
an enhancement of the ISW on very large scales (i.e. en- 
hanced decay of gravitational potentials) and on small 
scales there is a suppression of ISW effect (growth of po- 
tentials). Then, at late times z < zrs these corrections 
invert themselves and ISW is suppressed on large scales 
and small-scale potentials decay. 

Turning now attention to the auto-power correction, 
Eq. (|42|) . we see that the prefactor multiplying Pj^ is 
always positive, whereas the second term only switches 
sign when 1 - / < or 1 - 3/ > 0. However since 
Q m < 1, the first bracket will never switch sign and will 
vanish at early times. On the other hand the second 
bracket remains negative until Q m rs 0.16. Given current 
constraints on f2 m0 ~ 0.25 the second bracketed term is 
always negative and on multiplying by Pgg, we conclude 
that it too is always positive. 



D. Evolution of density power spectrum 

Fig. [3] shows the evolution of the nonlinear matter 
power spectrum measured from z = 1 to the present 
day. The power is plotted from the fundamental mode, 
k = 2tt/L « 0.005 /iMpc -1 to half of the Nyquist fre- 
quency of the mesh &ny = itNg/L ss 2ft,Mpc _1 , where 
we use N g = 1024 for all transforms. Above this fre- 
quency the power in the Fourier modes is affected by 
aliasing from smaller scales [55J. In the top panel we 



show the mean ensemble averaged absolute power from 
the simulations at each epoch, colored points with er- 
rors. On the largest scales k < 0.1, the power grows by 
a factor of ~ 2 from z = 1 to 0, and there appears to be 
very good agreement with the linear theory predictions 
on these scales (colored solid lines). On smaller scales 
the power is significantly amplified. 

In the bottom panel of the Fig. [3] we take the ratio of 
the data with the linear theory, and to see clearly the 
effects for each snapshot, we offset the curves by 0.1 in 
the vertical direction, with the solid colored lines being 
the baseline for each corresponding snapshot. We see 
that there is a small (» 2 — 3%) suppression of power 
at late times for modes 0.05 < k < 0.1 [ h Mpc -1 ], this 
is termed the 'previrialization feature' [54l l56l. |57|. On 
smaller scales (k > 0.1 /iMpc -1 ) the power is strongly 
amplified, compared to linear theory. In this panel we 
also present the predictions from the standard PT (de- 
scribed in jjlVCp and we see that it qualitatively captures 
the trends in the data. However, in closer detail, we see 
that the PT over estimates the power on smaller scales 
and that the predictions become progressively worse at 
higher expansion factors and higher wavenumbers. 



E. Evolution of the pseudo-momentum spectra 

In Figs S] and [5] we present the pseudo-momentum- 
density cross- and pseudo-momentum auto-spectra, re- 




0.1 
k[h/Mpc] 

FIG. 3: The time evolution of the nonlinear CDM den- 
sity power spectrum as a function of wavenumber. Top 
panel: colored points denote the absolute power and error 
bars are on the mean and are determined from the ensem- 
ble of simulations. The thin lines denote the linear theory 
and from top to bottom results are for expansion factors: 
a = {1.0, 0.93, 0.87, 0.76, 0.66, 0.62, 0.57, 0.54}. Bottom panel: 
the ratio of the power spectra with respect to the linear theory 
prediction. The thick solid lines denote the predictions from 
the nonlinear Eulerian PT. Note that for clarity the measure- 
ments have been offset by 0.1 in the vertical direction. 



spectively. Again the top part of each figure shows the 
absolute power and the bottom the ratio with respect to 
the linear model (c.f. Eqs [28] and HSJ . Note that the 
spectra are amplified relative to the density spectrum, 
and that on large scales this boost is well captured by 
multiplicative powers of af(a)H(a). In addition, we find 
that on very large scales, the momentum spectra also dis- 
play a previrialization feature and that the suppression of 
power appears to be deeper in both cases. Furthermore, 
on smaller scales the nonlinear amplification, which oc- 
curred at around k ~ 0.1/iMpc -1 for Pss, appears at 
larger scales in both cases, with the P uu , strongly am- 
plified by k > 0.7 h Mpc -1 . We compare these measure- 
ments with the predictions from standard PT and find 
for Pg^ reasonably good agreement on very large scales 
and an over prediction on smaller scales. However, for 
Puu t ne agreement is much better. 
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0.01 0.1 1 



k[h/Mpc] 

FIG. 4: The evolution of the pseudo-momentum-mass density 
cross-power spectra from z — 1 to 0, as a function of spatial 
wavenumber. Points and lines are as for Fig. [3] 




0.01 0.1 
k[h/Mpc] 



FIG. 5: The evolution of the pseudo-momentum auto-power 
spectra from z = 1 to 0, as a function of spatial wavenumber. 
Points and lines are as for Fig. [3] 



F. Evolution of the $ power spectra 



Having examined the individual components of the 
spectrum we may now sum them together with 
weights as given by Eq. (|23|) . Following Seljak [l5[ and 
Cai et al. [1711, we introduce the dimensionless and re- 



scaled form of Pj 



A| 4 (fc) 



4tt fc 3 P M (fc) 
~ (2tt) 3 [T{k)H{a)/a] 2 ' 
k 3 r„ 2P uS {k) 



2tt 2 



Pss(k) 



Puto (k) 



H(a)a(t) H 2 (a)a 2 (t) 



(44) 



(45) 



Fig. [5] shows the evolution of the ensemble averaged A| • , 
with errors on the mean. The top panel shows the ab- 
solute spectra for the 10 snapshots from z = 1 to the 
present day. Also shown as the thin solid lines are the 
predictions from the linear theory as given by Eq. (|30[) . 
Again, there appears to be good agreement on large 
scales, and nonlinear amplification on smaller scales. The 
bottom panel presents the ratio with respect to the lin- 
ear model, again we have offset different epochs by 0.1 in 
the vertical direction for clarity. There is clear evidence 
for nonlinear amplification of the spectrum on the very 



largest scales, and relative to linear theory this becomes 
increasingly more important at higher redshifts, as ex- 
pected. Indeed by k — 0.03 h Mpc -1 and at z ~ 1.0 the 
power is more than 10% in excess of the linear theory 
prediction, whereas at z = 0, a 10% amplification is only 
achieved by k ~ 0.07 h Mpc^ 1 . Here we also show the 
predictions from the NLO PT calculation from Eqs (j4"0|) 
and (|42|) , and we note a startlingly good agreement at all 
epochs. 

That the nonlinear effects become increasingly impor- 
tant at higher redshifts follows directly from the fact that 
1— /(a) — > as a <C 1. In this case, the only contribution 
to the spectrum comes from the nonlinear Rees-Sciama 
effect, and in the limit a — > it is given by 



&U(k)^—[l~2f(a)fP™(k) 



(46) 



On comparing our results with Fig. 1 from Cai et al. 
(T7I . we find qualitatively good agreement. However, on 
the largest scales their spectra do not appear to repro- 
duce the linear theory at high precision. The excess sig- 
nal that they find compared to the linear theory, we be- 
lieve, is a result of using the approximation / w fi„ 6 . 
Some of the discrepancy may also be due to cosmic vari- 
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d0/dt 



d0/dt 



FIG. 6: Evolution of the $ auto-power spectra from z = 1 to 
in dimensionless units, as a function of spatial wavenumber. 
Again points and lines are as presented in Fig. [3] 



ance, since they only show results for a single simulation. 
In that work the authors also proposed a nonlinear cor- 
rection formula for P^,^,, which has two free parameters. 
Since the PT has no free parameters, and as it provides 
an excellent match to the data for the scales k < 0.1 we 
consider our approach a sufficient description on these 
large scales. Such fitting would most likely be necessary 
on smaller scales for good agreement, but these scales 
are of diminishing importance for the calculation of the 
CMB Ci spectrum for I < 100. 



V. RESULTS: IMPACT ON CMB SPECTRUM 

The CMB temperature fluctuations arising from the 
ISW may be decomposed using a spherical harmonic ex- 
pansion, and the amplitude of each harmonic can be writ- 
ten as [llllll, 

d 3 k 



aj m = I —-Y£ m (k)A?(k, a ls ) , (47) 



with, 



2a 



At(k)= / d X -Mk X Mk,x(t)) , (48) 
Jo c 



0.10- 



10 

Multipole I 



FIG. 7: Relative error ([Cf 



Cf xact ]/Cp xact ) between 



the Limber approximation and the exact Ci computation of 
the CMB angular power spectrum. Results are shown for 
k — (I + 1/2)/ Da replacement. 



where we transformed Eq. (|TJ) to comoving geodesic dis- 
tance x an( i Xmax is the distance from the observer to 
the surface at which the ISW first becomes significant. 
The power in the harmonic multipoles may be calculated 
using the standard methods, and the ISW temperature 
spectrum may be written: 



CT T = - I dkk 



dxidX2ji(kxi)ji(kx2) 



4aia 2 



(49) 



In the limit (I > 10) we may use the Limber approx- 
imation to simplify the above integrals [see for example 
El, EE HH ■ Assuming that only modes transverse to the 
line of sight contribute to the signal and also that the 
power spectra are slowly varying functions of k, then the 
orthogonality of the spherical Bessel functions gives, 



dkk 2 ji(kxi)ji{kx2)'P a {k, Xi, X2) 



k 5 D (xi -X2) 



I 



Da(xi) 



(50) 



where Da(cl) is the comoving angular diameter distance 
(Da(ci) — x( a ) f° r fl & t space). On applying this approxi- 
mation the above expression reduces to the simple form: 



4V 
dx ~~&P i>& 



k 



W X )^ ;(51) 



./aP M (* = ^,o ; H[a)x 



1 



■(52) 
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1 10 100 



Harmonic multipole 1 

FIG. 8: Top panel: Angular power spectrum of CMB tem- 
perature fluctuations as a function of harmonic multipole. 
ISW contributions: solid green shaded region - linear the- 
ory with encompassing 1-a Gaussian error domain; red dash 
line - nonlinear PT; ensemble average TV-body measurement 
- blue triple dot dash curve. The CMB primary anisotropy 
spectrum is given by the magenta dot-dash curve and its 1-a 
Gaussian error domain by the solid magenta shaded region. 
Bottom panel: ratio of the nonlinear ISW spectra to the linear 
theory spectrum. While this can exceed the 1-sigma error of 
ISW alone (as shown), it is always less than the 1-sigma error 
band of the overall CMB (not shown), hence the nonlinear 
effects are not detectable in the CMB power spectrum. 



Fig. [7] compares the Limber approximate expressions 
for the angular power spectrum of temperature fluctua- 
tions with the exact spherical harmonic line-of-sight in- 
tegration. On scales I < 10, the Limber approximation 
is clearly poor with relative errors being > 10%. The 
transformation k = (I + 1/2)/ 'D a, as suggested by Ho 
et al. j27j, Loverde and Afshordi (5t|, improves the ap- 
proximation, but the errors still remain large. However 
for I > 10, the error is reduced and by / = 20 it is of 
the order ~ 5% (for a detailed discussion on the va- 
lidity of the Limber approximation for different power 
spectra, see Appendix ■ We shall nevertheless adopt 
the Limber approximation for our theoretical analysis, 
but note that if significant effects are apparent on multi- 
poles I < 30, then only a full spherical harmonic analysis 
will give robust results. However, this would necessarily 



involve computing the unequal time correlations of the 
Fourier modes of ^(k, t). 

Fig. [5J shows the results for the Limber approximated 
ISW temperature angular auto-power spectrum. We 
scale the Ci spectrum by 1(1 + 1) in the usual way and 
restrict our attention to angular modes / < 100. In the 
upper panel of the figure we compare three predictions: 
the linear theory calculation with 1-a cosmic variance er- 
rors, denoted by the green shaded region; the nonlinear 
PT, denoted by the red dash line; and the mean measure- 
ment from the A^-body simulations blue triple dot-dash 
curve. The 1-a green shaded error region was computed 
using the Gaussian variance formula: 



cr ~v/ sky 2*+i ' (53) 

where / s k y is the fraction of sky covered, and we shall take 
this to be of order unity. The estimate of the Ci spec- 
trum from the Af-body simulations was obtained by the 
following prescription: we first made an array of the mea- 
sured V^(k) spectra and divided this through by the lin- 
ear theory ISW power spectrum at that epoch. On very 
large scales the ratios all asymptotically approach unity 
and so the only evolution that remains to be modeled is 
the higher k domain. To do this, we employ the bi-cubic 
spline routine [55| and interpolate the spectra in log 10 [a] 
and log 10 [k] . Note that on scales greater than the funda- 
mental mode of the simulation cube the bi-cubic spline 
gives unity and we recover exactly linear theory. We em- 
phasize the importance of this step, since otherwise the 
Cj T spectra will be significantly reduced for I < 10, ow- 
ing to the finite volume of the simulations. Note that in 
order to avoid extrapolating the bicubic spline fits into 
regions where we have no measured data, the upper red- 
shift limit of the Limber integrals was set to z = 1. We 
have tested that this does not change our results in any 
significant way, by computing the PT out to z = 5. 

In Fig. [5J we see that all three theoretical predictions 
converge for / < 30, however for I > 30 we find enhance- 
ment of the signal for both the PT and A^-body results 
and that these agree to high precision, in agreement with 
expectations from Fig. O By I = 50 they both show be- 
tween ~ 10—15% increase in the power. We also show the 
CMB primary anisotropy power spectrum as the black 
dot-dash line, with the magenta shaded band giving the 
cosmic variance errors, Eq. (|53[) . The primary Ci spec- 
trum was obtained using the cmbf ast routine with cos- 
mological parameters to match those of the zHORIZON 
simulations. Note that, by default, this spectrum already 
includes the linear ISW effect. 

Comparing the primary with the ISW signal, we see 
that at I — 30 the primary signal is two orders of magni- 
tude larger, and so the nonlinear enhancement at these 
multipoles will induce changes to the CMB spectrum that 
are 1%. While the nonlinear effect exceeds cosmic 
variance in ISW for I > 50, it never exceeds the cosmic 
variance from the total CMB, since ISW contribution 
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to CMB decreases with I. Our findings are consistent 
with earlier results [12, E3] > but are established with im- 
proved precision. We therefore do not expect large-scale 
nonlinear evolution of the gravitational potentials to be 
responsible for any anomalies in WMAP angular power 
spectrum. 



VI. ISW-DARK MATTER 
CROSS-CORRELATION SPECTRUM 

Having discussed the ISW auto-correlation spectrum 
we now move on to discussing ISW correlation with the 
density field. We begin with the dark matter density 
Cf s . This can be observed by cross-correla ting the CMB 
with the weak lensing signal of galaxies [6Ct |6l| , the weak 
lensing of 21cm transitions [62| or the weak lensing of 
CMB itself (with information encoded in CMB bispec- 
trum) [33|, [63[ . In addition, there are a number of ad- 
vantages to be gained from studying this: firstly, there 
exists an "alternative" method for estimating P$ s , and 
this provides us with an independent check on our "stan- 
dard method", described in ffTTJ secondly, owing to the 
larger number of dark matter particles, the effects of shot 
noise on the spectra can be better assessed, and as we will 
show for the alternate method, more easily corrected for. 



A. Alternative estimator for P^, s 

Our alternative approach to estimating P^ s can be un- 
derstood as follows: Consider the ensemble average of 
the product of 5(k) and using Poisson's equation 

we may rewrite this as, 



P 5i (k,a) = ^(<5(k, a )<f*(k,a) 



-a[^(A:)]- 1 ^($(k ) *)$*(k ) a) 



= -a[^(fc)]- 1 P (b$ (k, 



(54) 



We now take advantage of the useful property that 

19 Id 
f*i(k,a) = -— P$$(k, a) = -aiJ(a)— P$$(k, a) . 

(55) 

Through further use of Poisson's equation, the last 
term in the above equation may be rewritten in terms 
of the density power spectrum, i.e. P$$(k, t) — 
[J-(k)/a] Pss{k, a). Putting this together, we arrive at 
the result [33j . 

P is (k,a) = --a H{a)T{k) — 

1 (9 In y 



-aH{a)y{k,a)T{k) 
2 a In a 



(56) 



where y = [Pss (k, a) /a 2 ] . This simple expression in- 
forms us that the ISW cross-correlation can also be es- 
timated from just two things: the matter power spec- 
trum and its evolution with time. We may check that 



the above expression is consistent with our previous re- 
sult (c.f. Eq. l2"4"|) . On assuming linear theory Pss{k, a) — 
D 2 (a)Pun(k), then we find 



dirty 
d In a 



2 if (a) - 1] 



(57) 



and on insertion of the above expression into Eq. ([52 
we recover our earlier result. 

The practical implementation of the above algorithm 
requires us to estimate the time derivative of the power 
spectrum, and we do this using the usual time-centered 
difference scheme: 



■Vi-i 



dh\y ^ 1 Hi. 

d log a yi [log a l+i - log a { 



(58) 



where yi = yifli) is the estimate at epoch o^. Note that 
since we employ a time-centered difference scheme, we do 
not show results for for z = or z = 1, the first and last 
epochs considered. 



B. Results: evolution of P^ s 

Fig. O compares the results for P^ s obtained from our 
standard method (c.f. Eq l24j) of solving the continuity 
equation (black points with errors), with our alternative 
method (colored points with errors). As was done for 
Pj>$ we have introduced a dimensionless and scaled form 
of the cross-power spectrum (c.f. Eq. I45|) : 



A 2 , 

(5* 



(fc): 



-in 



(2tt) 3 [F{k)H(a)/a] 



_*1 

2^ 2 



Pss(k) 



Ps^k) 



H{a)a{t) 



(59) 
(60) 



The left panel of Fig. [S] shows that the two independent 
approaches produce results that agree to high precision. 
We are therefore confident that both methods are consis- 
tent and implemented correctly. 

The top panel compares the spectra estimated from 
the simulations (points with error bars) and the linear 
theory predictions (solid lines). The lower panels show 
the ratio with respect to the linear predictions. There 
are a number of important features that we draw atten- 
tion to: firstly, rather than nonlinear effects becoming 
increasingly prominent with time, we see that they are 
stronger at earlier times and on larger scales. The expla- 
nation follows our earlier discussion of Fig. [51 and owes to 
the fact that the linear ISW effect switches off as a — > 
and 1 — f(a) — > 0, leaving only the RS contributions 

(c.f. gEE). 

Next, we note that there is a sign change in the spectra 
as one goes from low to high k. Since we plot the absolute 
value of the power the sign change is understood to be 
the point where the signal drops to zero and bounces 
back up. The scale at which this sign change occurs is a 
function of time, and it appears on larger scales at higher 
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FIG. 9: Comparison of P^ s estimates obtained from standard method (continuity equation) and the alternate method (time 
derivative of Pss, c.f. Eq.l56[l. Colored points with errors denote the standard method, and black points with errors denote the 
alternate method. Left panel: no shot noise correction. Right panel: shot noise correction applied to the alternative method. Top 
sections show absolute power, lower sections show the ratio with respect to the linear theory, and for clarity the results at each 
epoch have been offset from each other by 0.1 in the vertical direction. The lines in the bottom panels denote the predicitons 
form the standard PT and from top to bottom results are for expansion factors: a — {1.0, 0.93, 0.87, 0.76, 0.66, 0.62, 0.57, 0.54}. 



redshifts [see also [TtJ • The sign change is due to the fact 
that the $ signal becomes dominated by the RS and BG 
effects. However, for the spectra with z < 0.3 we see no 
sign change over the fc-range that we consider. Moreover, 
unlike the lower redshift epochs we see an amplification 
relative to the linear theory. This means that, at late 
times in LCDM model, nonlinear evolution can actually 
enhance the decay of gravitational potentials, consistent 
with our earlier discussion of the PT (c.f. ^IV CI) . Further 
support for the PT interpretation of this phenomenon 
comes from the fact that if one considers the results at 
high redshift, then around k ~ 0.1/iMpc -1 there is a 
small amplification of power with respect to the linear 
model. 

In the discussion so far we have neglected the issue of 
the discreteness correction due to finite number of dark 
matter particles. It is unclear how to apply the shot noise 
correction to the momentum density. However, since we 
know the shot noise correction on the dark matter power 
spectrum is Pss — * Pss ~ where n is the number 
density of dark matter particles, the discreteness effects 
may be accounted for more easily when using Eq. (|56[) . 
Fig. [HI right panel, shows the results obtained from this 



procedure. Whilst we see that the correction reduces the 
spectra by a small amount for all k > 0.2, we neverthe- 
less see that both the small-scale late-time amplification 
and early-time large-scale amplification of the P^ s re- 
main significant. We are therefore led to conclude that, 
nonlinear evolution may lead to a small enhancement of 
the ISW in the LCDM model. 

Comparing our results with the measurements of P s ^ 
from Cai et al. fl7j |. we observe that these authors find 
no such late time amplification. Owing to the fact that 
we have provided two independent methods to obtain 
the estimates, and since we have a significantly larger 
total simulation volume (~ 12 times larger) furnishing 
smaller errors, we believe that our result is robust. In 
the next section we shall investigate whether selecting 
highly biased regions may influence these results further. 



VII. ISW-HALO CROSS-CORRELATION 
SPECTRUM 

The cross-correlation between $ and a density tracer 
field is more easily observable than with the density field 
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itself, which so far is limited because of the small area 
or large errors in the weak lensing reconstruction. One 
consequence of this is the added complication of need- 
ing to understand the bias relation - the mapping from 
the tracer population to the underlying dark matter den- 
sity. In this Section we shall explore whether the cross- 
correlation of $ with cluster- and group-scale dark matter 
haloes, measured in the zHORIZON simulations, between 
z = [0.0, 1.0], changes the ISW signal in any significant 
way, beyond a linear bias. From the assumption that 
all galaxies reside in dark matter haloes, it follows that 
the large scale clustering properties of any galaxy sample 
are a weighted average of the halo clustering statistics. 
Consequently, studying the halo-ISW cross-correlations 
should provide representative results for a number of 
plausible surveys. In particular, while we are limited by 
the mass resolution of our simulations so that our analy- 
sis only applies to biased halos with bias b > 2, we note 
that most of the data sets used for ISW detection so far 
are based on strongly biased tracers Giannantonio et al. 
[26| . Ho et al. (27[, so our results are applicable to these. 



A. Linear Theory 

In nearly all ISW studies to date the bias has been 
assumed to be not only constant in space, but also in 
time. As discussed in Ho et al. 27J and more recently 
Schaefer et al. [64], if one wishes to go beyond detection 
and constrain cosmological models with the ISW, then 
it is likely that this over simplification will introduce a 
bias in the recovered parameters, especially when redshift 
selection functions are broad. The next simplest scenario 
is a time-dependent linear relationship: 



5 Q (x, a) — &"(a)(5(x, a) 



(61) 



where S a — > {g, h, c, . . . } denotes the tracer type, e.g. 
galaxies, haloes, clusters etc., bf (a) is a linear bias pa- 
rameter that varies in time but is independent of scale. 
In this case the ISW cross-spectra and biased tracer auto- 
spectra may be easily computed as (cf. Eq. I31[) : 



H(a) 



(1-/W) 



Pjf(M) £62) 

(63) 



B. Nonlinear theory for the bias 

Several recent theoretical and numerical studies of the 
bias of dark matter haloes [13, [57!, 65], have revealed that 
the linear model is only likely correct on asymptotically 
large scales. These predictions have been confirmed by 
several observational studies of the relative bias of differ- 
ent galaxy populations [H, [13, [HI ■ In Smith et al. [EU 
it was shown that the scale dependence of halo bias was 
a strong function of scale for k > 0.07 /iMpc -1 . In that 



work a physically motivated analytic framework was de- 
veloped to model these scale changes. A similar approach 
was independently developed by (M 70t. The model uti- 



lizes a nonlinear local bias model 
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8 a (x,a) = ^2 



F(x,a)-(5"(x))] 



(64) 



where the constant term from the Taylor expansion was 
rewritten as 6g = — ^,- =1 bj(a) (£ J ) The density 
field may be expanded using the PT series expansions 
from HIV CI As was shown by Smith et al. 54|, if one 
transforms to Fourier space and collects terms to a fixed 
order, then the density field of the biased tracers may be 
written as a fluctuation series of the form: 



6 a (k,a\R) 
[S a Qt\a,R)] n 



J2[D(a)] n [S a (k,a\R)] n ; (65) 



(2tt) 3 ™- 3 
xF£(<ii,...,q n \a,R) , 



(66) 

where [<5 Q (k|a, R)] n is the nth order perturbation 
to the biased tracer density field. The functions 
F"(qi, q n |a, R) are the bias tracer PT kernels, sym- 
metrized in all of their arguments. The kernels are de- 
scribed in Smith et al. [54[ . Thus equations (|6"5"|) and (|6"6"1) 
can be used to describe the mildly non-linear evolution 
of the biased fields to arbitrary order in the dark matter 
perturbation. There is a subtlety that we have skipped 
over: in order to facilitate the Taylor expansion of the 
biased field it was necessary to filter on a scale R, and 
hence all of the kernels depend on the filter scale. To 
remove the filter dependence we adopt the renormaliza- 
tion scheme suggested by McDonald [i^, [7(| • The down 
side of this, is that the parameters may not be derived ab 
initio, but must be obtained through fitting to measured 
data and we shall do this in the following section. 

Using these relations, along with McDonald's renor- 
malizations, we find that the ISW-biased density tracer 
cross- and auto-power spectra may be written: 



P 



NL 



(k,a) = P^(k,a)+Pll oop (k,a) 
P^(k,a) = P^(k,a)+P^(k,a) 
where the loop corrections are given by, 



(67) 
(68) 



qILoop 
a<i> 



H{a) 



HQ — [1 - V{a)] P£7(fc, a) ; (69) 



P 



lLoop 



b%iPll°° v + b a Ra A{k,a) 



(70) 



\h a l 2 
°R,2J 



Pa L a ° op = 2bl 1 b^ 2 A(k,a) + ^f-B(k,a)+N%m) 
and where we have introduced the auxiliary functions: 

A(k,a) = y^P L i„(9)PLin(|k-q|)F 2 (q,k-c072) 
B(k, a) = J ^3 P Li n(q) [P L in(|k - q|) - P L in(<7)l73) 
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TABLE II: Halo mass classes, and number densities. 

Mass Range n(z — 0) n(z — 1) 
[xHF/T^Mq] [h- 1 Mpc\ 3 [/t^Mpc] 3 



Bin 1 
Bin 2 



[1.50 < M < 10.0] 3.5 x 10" 4 1.8 x 10 
[10.0 < M] 2.5 x 10~ 5 3.3 x 10 



In the above equations we introduced the renormalized 
bias parameters b^^a) and the renormalized constant 
power term iV^(a). This may be thought of as an arbi- 
trary white noise contribution. 

Before moving on, we notice that the sign reversal 
property of the nonlinear cross-power spectrum of $ with 
mass density, remains unchanged. This owes to the fact 
that bfi t 2 changes the scale at which the loop corrections 
transit from large-scale power suppression to small-scale 
enhancement (provided loop corrections are small com- 
pared to linear theory). 



C. Renormalized halo bias parameters 

In order to use the nonlinear bias model we require the 
time evolution of the bias parameters, b\ and b 2 . These 
can be estimated directly from the simulations in the 
following way. Firstly, we divided the haloes at each 
expansion factor into two classes: (Bin 1) group scale 
dark matter haloes and (Bin 2) cluster scale dark matter 
haloes (see Table ILT1 for details). 

These mass bins can be faithfully traced within our 
simulations out to z = 1. We choose fixed mass bins 
at all epochs for simplicity, but a reasonable association 
can be made between these halo bins and tracer popula- 
tions such as Luminous Red Galaxies (LRGs) or clusters. 
Then for each realization we compute the power spectra: 
^hh, PhS, P\lu), and -P h $, for all of the snapshots from 
z = [0,1]. The renormalized halo bias parameters were 
then directly estimated from the P^s data in the follow- 
ing fashion. Firstly, we fit for b\ l on the largest scales 
using an inverse variance estimator of the form: 



J R,l 



(74) 
(75) 



where b\(ki) = PhS,i/Pss,i and with [ki,k 2 ] = 
[0.0, 0.05] /iMpc -1 . Note that we assume that there is 
little covariances between k bins on these large scales. 
Having obtained b^ x , we next obtain our estimate for 
b\ 2 ■ Our estimator has exactly the same form as the 
above equations except for the fact that b\ 1 

a bl -> a b2 and that [k u k 2 ] = [0.05, 0.2] h Mpc -1 
important quantity to specify is, 



b\ 2 and 
The 



b h 2 {h) = 



(76) 




FIG. 10: Renormalized bias parameters bn,i and 6h,2 as a 
function of redshift. Red (thin) and blue (thick) lines and 
points denote results from Bin 1 and Bin 2, respectively. The 
symbols are: estimates - solid points; bn,2 estimates 

- stars. The dotted lines show the bias evolution model of 
Eq. |78]l. 



It will also be useful later on for us to predict Phh, 
and to do this we are required to additionally estimate 
the renormalized shot noise term: JV^(a). This may be 
obtained directly from our estimates of Pth,i along with 
Eq. (USD and by using Eq. $7^, but with b^ R1 -> 

t&i -> and with [fci,fc 2 ] = [0.0, 0.2] /iMpc -1 
estimate per mode is 



N\ and 
Our 



p 



hh,i 



[b h n,i?Pss(k) 



[b\ 2 fB{ki 



(77) 



Fig. [TU] shows the time evolution of the best fit renor- 
malized halo bias parameters. As is evident from the fig- 
ure, the values of b\ 1 for the two samples decrease with 
increasing time. This is qualitatively consistent with the 
halo bias evolution that emerges from Extended Press- 
Schechter formalism and the Peak-background split ar- 
gument (dotted lines), where linear halo bias decays with 
time as [zl, Lz3, Elj : 



[61(a) - 1] = D(a )/D(a) [h{a ) - 1] 



(78) 



However as the figure clearly shows the actual measured 
halo bias evolves much more strongly as a function of 
redshift. We also note that the values of b 1 ^ 2 are also 
similarly consistent with this theory, which predicts that 
fe'p 2 < for haloes around M* (t) (the characteristic non- 
linear halo mass at that epoch a(M* , t) = 1), and that 
b\ 2 > for haloes with M > M*(t) [76j . 



D. Results: Evolution of halo— density spectra 

Fig. [11] shows the evolution of P^s in the simula- 
tions. The left panel presents the results for haloes in 
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k[h/Mpc] k[h/Mpc] 

FIG. 11: Evolution of the h—8 cross-power spectra as a function of spatial wavenumber from z = 1 to 0. Left panel is for haloes 
in Bin 1 and the right panel shows results for haloes in Bin 2. Points and lines are as presented in Fig. [3] 




0.01 0.1 1 0.01 0.1 1 



k[h/Mpc] k[h/Mpc] 



FIG. 12: Evolution of the h-cu cross-power spectra as a function of spatial wavenumber from z = 1 to 0. Left panel is for haloes 
in Bin 1 and the right panel shows results for haloes in Bin 2. Points and lines are as presented in Fig. [3] 
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Bin 1 and the right Bin 2. The top sections show the 
absolute power and the lower sections show the ratio 
with respect to the linear theory predictions. We ob- 
serve that all spectra exhibit a strong scale dependence 
relative to the linear theory and that the departure is 
characterized by a suppression of power on large scales 
(k > 0.05 /iMpc -1 ) followed by power amplification on 
smaller scales (k > 0.1 /i Mpc -1 ), and this is exhibited 
in both mass bins and at all times. The highest mass 
bin exhibits the strongest amplification with scale, by 
k ~ 0.1 ft. Mpc -1 the spectra are 10% in excess of the 
linear theory, whereas Bin 1 shows a slightly stronger 
large-scale power suppression. In the lower sections of 
each panel in Fig. [11] we also show the predictions of the 
nonlinear renormalized bias model from § I VII Bl and we 
find surprisingly good agreement over all of the scales of 
interest. We note that for Bin 1 on the smallest scales 
k > 0.5 h Mpc -1 , the predictions appear to drop dra- 
matically to zero. However, for the computation of the 
Ci we expect that this theoretical accuracy will be suf- 
ficient. This owes to the fact that the spectra shown in 
Fig. [11] will all be premultiplied by J-(k) oc fc -2 and so 
will be suppressed relative to larger scales. We note that 
the scale dependence of the halo cross-power spectra were 
investigated by Smith et al. and we confirm the basic 
results presented in that study. 



E. Results: Evolution of halo— momentum spectra 

Fig. rT3] shows the evolution of the halo-pseudo- 
momentum cross-power spectra as a function of scale. 
Again left and right panels are for Bins 1 and 2, respec- 
tively. As expected from our investigation of P$u, we 
again see nonlinear features in these spectra, and that 
they are more enhanced relative to those in the P^s spec- 
tra. This can be inferred through considering the ratios 
of the spectra with respect to linear theory (bottom sec- 
tion of each panel). In particular, we note that for Bin 1 
and the z — snapshot, the large-scale suppression fea- 
ture is of the order ~ 5% at k ~ 0.7 h Mpc -1 , in contrast 
to ~ 2% supression in P^s- Again, in the lower pan- 
els we over plot the predictions from the renormalized 
bias model and the agreement is again good, although 
for k > 0.2 h Mpc -1 small deviations of the model from 
the data are more apparent. Also, the predictions for Bin 
1 drop to zero at higher k, and this occurs for the reasons 
previously noted. 



F. Results: Evolution of halo— $ spectra 



In Fig. [13] we combine the power spectra from the pre- 
vious two subsections to explore the evolution of P,i . As 
was done for the analysis of P&& and P s ^ we introduce a 
dimensionless and scaled form of the biased cross-power 



spectrum (c.f. Eg. [45]) : 

k 3 r 



h<t> 



2tt 2 



H{a)a{t) 



(79) 



The top panels compare the spectra estimated from the 
simulations (points with error bars) and the linear the- 
ory predictions (solid lines). The lower panels show the 
ratio with respect to the linear predictions, and the lines 
show the predictions from the renormalized nonlinear 
bias model. As was the case for our investigation of A \ . 
(c.f. IVIBj) . departures from linear theory are increasingly 
apparent as one considers higher redshifts. In addition, 
there is a sign change in the spectra as one goes from 
low to high fc. The explanation again follows our ear- 
lier discussions surrounding Figs [6] and [9] On comparing 
these results for the haloes with those for the dark mat- 
ter, Fig. [5] we find that the scale at which the spectra 
switch sign becomes larger with increasing bias. 

Considering the small-scale, late-time ISW boost rela- 
tive to linear theory, we see that for the haloes at z = 
the signal is stronger as bias increases. However, we also 
note that the amplification is present for the Bin 1 halo 
sample by z ~ 0.3, compared to the Bin 2 sample where 
it is absent by z > 0.1. This result means that, at late 
times in LCDM model, nonlinear evolution can enhance 
the decay of gravitational potentials and that the rate of 
decay also depends on the environment. Again, this re- 
sult naturally emerges from the PT (c.f. fllVCp . although 
as is shown in the figure, the PT struggles to capture the 
measured spectra precisely. In the next section we shall 
investigate whether these nonlinear effects are sufficiently 
large to impact the ISW-density tracer C/'s. 



VIII. CMB-LSS ANGULAR POWER 
SPECTRUM 

A. Theory 

We now turn to the calculation of the ISW-biased den- 
sity tracer angular power spectrum. As described in |jV] 
for the ISW auto-spectrum, we may also decompose the 
projected fluctuations in our biased density tracer into 
spherical harmonics. To do this, we define the 2D biased 
density field as the weighted projection of the 3D density 
field along the line of sight and in a cone of solid angle 
d£l. This we may write as, 

ST(0) = r d X D\{a)q(x)5T{D A (x)lx) , (80) 

where q(x) is a radial weight function, which is normal- 
ized such that 



d X ^D 2 A ( X )q(x) = 1 • 



(81) 



To proceed we must specify q(x)- For a typical mag- 
nitude limited survey the weight function would be 
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FIG. 13: Evolution of the dimensionless scaled h-<I? cross-power spectra (A^) as a function of wavenumber from z = 1 to 0. 
Left and right panels show results for haloes in Bin 1 and Bin 2, respectively. Points and lines are as presented in Fig. 



l(x) — n i> L x )/Ntot where n(> L x ,x) is the space 
density of galaxies above the flux limit at a given red- 
shift, and -/Vtot is the total number of galaxies, and 
so D\(x)q(x) tx dN(z)/dz the number redshift distri- 
bution. Therefore, in turn, one is required to specify 
a model for the redshift distribution [see for example 

mm. 

Since we are more interested in precisely quantifying 
the importance of nonlinear contributions to the cross- 
correlation signal for biased tracers, which we can mea- 
sure directly at all epochs in the simulations, we shall 
therefore forgo attempting to fabricate certain aspects 
of a real galaxy survey - this level of detail may con- 
fuse interpretation. Instead we shall take a more simpli- 
fied approach: we assume that, above some fixed mass 
threshold, there is one and only one galaxy (perhaps an 
LRG) per dark matter halo; that the mass threshold is 
independent of redshift; and that we may construct a vol- 
ume limited sample of these objects from z = out to 
z = 1. This last condition implies that there is a tight 
relation between the mass threshold of the host halo and 
the luminosity threshold for the carefully selected target 
galaxy. Our model galaxy survey is therefore equivalent 
to a target sample of haloes above some fixed mass from 
redshift z — to 1. Hence, we shall write the weight func- 
tion, q(x) = n(> M,x)/N T oT(Xi,Xj), where n(> M,x) 



is the cumulative number density of dark matter haloes 
with masses greater than M at time i(x); and where by 
our normalization condition, for a redshift shell between 
Zi and Zj we have 

N T OT(Xi,Xj)= d X ^D 2 A ( X ) dMn(M; X ) ■ 
Jxi Jm 

(82) 

In the above, n(M,x), is the differential halo mass func- 
tion at time t(x) and Xi = x( z i)- Figure [bH shows the 
redshift distributions of our mock target samples, in the 
two mass bins and as a function of redshift. In the figure 
we have introduced the new weight function 

U tJ (x)^4nD 2 A ( X )e tJ (x) J™ dM n{M ' x) , (83) 

Jm N T oT[Xi,Xj) 

where Q lJ {x) = [©(x - Xi) - ®(x - Xj)], is the top-hat 
function with being the Heaviside step function. 

The multipole amplitudes of the biased density tracers 
may therefore be written, 

«L = H)'4tt J ^Y* m (k)A^k, Xl ,Xj) , (84) 
with, 

Afi^xuxj) = r rfxn y (x)jKfcx)^h D (fc,x) • (85) 



20 




and for the halo auto-power spectrum we have 



FIG. 14: Mock LRG/cluster normalized number redshift dis- 
tributions as a function of redshift. Top panel shows results 
for intermediate mass host haloes (Bin 1), and lower panel 
results for cluster mass host haloes (Bin 2). Note that here 
we show the normalized distributions over the entire sample 
range z = 0-1. The blue vertical dash lines show the 5 red- 
shift bands for which we compute the cross-correlations, and 
note that we renormalize the distribution for each band. 



Following Eq. (|4"9")l , the cross-angular-power of the ISW 
temperature fluctuations and the projected density trac- 
ers may then be written: 



9 r rx* 
C f T = ± I dkk 2 

7T 



dxidX23i{kxi)ji{kx2) 



x -^T U t](x)P h i(k;Xi,X2) 



(86) 



Under the Limber approximation (c.f. this expres- 
sion reduces to: 



C, 



hT 



i S a{Xm .J lna ^ {a)p ^{ k ^ D^Y a . 

x — 1 - ; ; (87) 



dlnaL T ? (a)Phh k = - . r ,i 

a( Xmax ) V D A{a) 

1 



H{a)x 2 {a) 



H(a) X 2 {a) 

In Appendix lAl we present a short investigation of the 
validity of the Limber approximation for predicting the 
ISW-LSS cross-power spectrum. There we show that the 
relative error is < 10% for I ~ 10 and that for I > 10 it is 
< 2%, and for a wide range of survey window functions. 
These results are consistent with the findings of Rassat 
et al. [25j for the 2MASS survey. Since we are inter- 
ested in scales I > 10, we shall therefore use the Limber 
approximated expressions. 



B. Results: ISW— biased tracer angular spectrum 

Figure [15] presents the results for the angular cross 
power spectrum for the ISW and haloes in Bin 1 (left 
panel) and haloes in Bin 2 (right panel). In each case 
we show the results for 5 narrow bins in redshift, and 
where for each bin we weight by the redshift distribu- 
tions presented in Fig. [JJJ The solid green lines in the 
figure denote the linear bias predictions; the red dashed 
lines correspond to our predictions from the nonlinear 
renormalized PT, as described in WII Bl and WII Ct and 
the blue triple-dot dash lines correspond to our bi-cubic 
spline fit to the ensemble average measurements of P h ^ 
from the simulations, and scaled by linear theory. 

In the figure we see that for both Bins 1 and 2 the peak 
of the angular power spectrum moves to the right and 
upwards as the mean redshift of the sample increases. 
The rightward shift is due to the fact that for a given 
physical scale the angular size decreases with distance, 
in this case the scale is the peak of the P h< j spectrum. 
The upward shift is more complex, if we were considering 
unbiased tracers then we would expect that the signal 
would drop with increasing redshift, owing to the fact 
that the ISW signal switches off and also the amplitude of 
the power spectrum is decreasing with D 2 (z). However, 
for a fixed mass range, the bias of the sample increases 
with increasing redshift (c.f. Fig. ITU]) . For the two host 
halo mass bins that we consider the bias evolves by a 
factor of ~ 2 from z = 0-1. 

Regarding the impact of nonlinearity on the predic- 
tions, we find that for I < 100 these are small, being at 
most < 10%. For Bin 1, the deviations are characterized 
by a several percent boost around I = 50, followed by 
a several percent suppression by I — 100. Whereas for 
Bin 2, the deviations are represented as a few percent 
suppression. For I > 100 the deviations are, in all but 
one case, characterized by a much more significant sup- 
pression, and the signal rapidly drops to zero. The case 
which does not conform to this picture is the lowest red- 
shift slice for Bin 1, here the signal estimated from the 
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FIG. 15: Angular cross-power spectrum of ISW effect and haloes as a function of spherical harmonic multipole I. Left panel: 
results for group scale dark matter haloes. Right panel: results for the most massive clusters. In each panel we show results 
for 5 equally spaced bins in redshift over the range: z — [0.0, 1.0]. The predictions are differentiated by line thickness: thick 
lines - low redshift; thin lines - high. The line styles denote: linear theory - solid green line; nonlinear PT - red dash line; 
bi-cubic spline fit to the simulation data - blue triple-dot dash line. Top sections of each panel give the absolute power, and 
the lower sections the ratio with respect to linear theory. The shaded regions represents the 1-a error domains per multipole 
of the linear cross-spectra, where the central redshifts z g {0.1,0.3,0.5,0.7,0.9}, correspond to the colours (red, green, blue, 
cyan, magenta). 



simulations appears to be boosted by ~ 10% at I ~ 100. 
Unfortunately, this amplification is not mirrored in the 
predictions from the PT, as also seen in Fig (TT5|) for the 
last four spectra. 

In Fig. [15] we also show the expected 1-a error do- 
mains (shaded regions) of the cross-spectra, computed 
from using the simple variance formula: 



A 2 (Q Th ) 



1 



1 



/sky (21+1) 



(89) 



As in the case for Cf T , we again find that the cosmic 
and sample variance errors dominate over the modeling 
errors on scales / < 100. 



C. Calculation of the S/M for biased tracers 

The result from Hernandez-Monteagudo [29[ is that for 
the ISW-dark matter cross-correlation up to 90% of the 
Signal-to-Noise (S/M) for the ISW comes from harmonic 
modes Z < 50. Here we shall assess whether sampling bi- 
ased density tracers can change these conclusions. From 



the last equation we write the S/M for the ISW-dark 
matter cross-correlation, at a given multipole I, as 



(S/M)* = / sky (2/ + l) 



T8\ 



(90) 



Similarly, this equation can be written for the halo dis- 
tribution, 



(S/M)* = Uy(2l + 1) 



L C TT C hh 



(91) 



In the above, no shot noise subtraction on the halo auto- 
power spectrum is assumed. We can define the cumula- 
tive S/M below a given multipole I as 



(S/M) [< I] 



EM- 



(92) 



This addition is legitimate only under full sky coverage 
(/sky = 1), since we assume that different multipoles are 
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FIG. 16: Top and bottom panels, S/N results for Bins 1 and 
2, respectively. Left panels: Squared S/N for each multipole 
I of the cross correlation of the CMB temperature with the 
most massive (Bin 2) halo population (solid lines) and the 
total matter density field (dashed lines). Results for differ- 
ent snapshots centered at z = 0.1, 0.3, 0.5, 0.7 and 0.9 are 
displayed in black, blue, green, red and orange colors, re- 
spectively. Line thickness decreases with increasing redshift. 
Right panels: cumulative S/N. Lines are as in left panel. 



independent. In the left panel of Fig.[16]we show (S/N)f, 
for the cluster-mass halo population (Bin 2, solid lines) 
and the matter density field (dashed lines) for the 5 dif- 
ferent redshift shells. The right panel of the figure shows 
the corresponding cumulative S/N below each multipole 
I. We note that the (S/N)f is flat for low multipoles, and 
declines rapidly with increasing I. The scale at which the 
turn down occurs is a function of redshift. For z = 0.7 
the turn down occurs at I < 20, whereas for z = 0.1 it has 
dropped by I ~ 10. From studying the cumulative S/N, 
we find that roughly 50% of the total S/M is achieved 
by I ~ 10, and that - 90% is achieved by I » 40 [c.f.H|. 
On comparing these results with the corresponding ones 
for the matter field (dashed lines), we find slightly lower 
values for the haloes. This may be atributed to the ad- 
ditional Poisson noise. Note that the redshift shell that 
gives the highest S/N is located at z = 0.7, and that the 
total S/N for it is of order ~ 7. 

Based on these findings, we conclude that it is highly 
unlikely that nonlinear evolution of the mass distribu- 
tion or nonlinearities in the scale dependence of bias can 



significantly affect the detectability of the ISW. 



D. Results: Cross-correlation coefficient 

Finally, we investigate the cross-correlation coefficient 
of the CMB temperature fluctuations and the halo sam- 
ples. The cross-correlation coefficient of two fields A and 
B is defined as, 



AB 



(jAAfjBB 



(93) 



Under the assumption of time independent linear bias 
we would have r Th — > r TS . Thus rf h does not depend 
on the bias of the tracer sample, nor the amplitude of 
the primordial power spectrum. Instead it provides di- 
rect information on the Dark Energy parameters and the 
curvature density: {£Ide, wq, fife, }• This approach was 
developed by Giannantonio et al. [26[ to obtain cosmo- 
logical parameter constraints from current CMB and LSS 
data [see also [77l . for an alternate method for removing 
bias, that uses CMB lensing.]. 

The validity of this analysis hinges on the fact that 
bag cancels out. However, since the bias is in fact time- 



dependent, we can only have 



„T6 



Adding to this 



the fact that the bias is scale-dependent it appears that 
such an approximation is unlikely to be robust, and espe- 
cially for LSS surveys with broad selection functions. We 
may test their conjecture by estimating r Th for several 
samples of biased tracers, and if we do not find that they 
match r TS within the same redshift shell, then the mod- 
eling should be deemed to be insecure. In that case one 
must include the redshift evolution of the bias, as done 
by Ho et al. [27| in their analysis of the NVSS sample. 

In Fig. [T7] we present the measured cross-correlation 
coefficients for the Bin 1 (left panel) and Bin 2 (right 
panel) halo samples and for the 5 redshift bins previously 
considered. The linear theory predictions are represented 
by the solid green lines and note that for these we use the 
time-dependent linear bias estimated directly from the 
simulations. In the figure we also present two different es- 
timates for the full nonlinear cross-correlation coefficient, 
estimated from bicubic spline fits to the measured spec- 
tra: the blue triple-dot dash curves show the results for 
the case where no shot noise subtraction was performed 
on the C; hh data; the red dash curves show the same but 
with the shot noise subtracted. We also show the dark 
matter-CMB cross-correlation coefficient, r TS , measured 
in the same redshift bins as for the haloes (magenta dot- 
dash curves). For the dark matter estimates, we used the 
selection function Hg M ( X ) = D\{ X )/ S% *xD\{x) ■ 

For these narrow redshift shells, Az — 0.2, we find that 
for linear theory, neglecting the evolution of the bias does 
not lead to significant errors. This can be seen from the 
middle panels of the figures, where we plot r Th /r TS (solid 
green line for linear theory) . However, for the actual mea- 
sured nonlinear r Th , we find that the scale-dependence 
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FIG. 17: Cross-correlation coefficient between CMB and haloes (r ) as a function of harmonic multipole Left panel 
shows results for group-scale haloes and the right for cluster-scale. Upper panels: results for 5 redshift bins centred on 
z = {0.1, 0.3, 0.5, 0.7, 0.9} of thickness Az = 0.2. Thick to thin lines denote low to high redshift halo samples. Line styles are: 
green solid line - linear theory; blue triple-dot dash line - bi-cubic spline fit to simulation data without shot-noise subtraction; 
red dash line the same but with shot-noise subtraction. The magenta dot-dash curve shows the result for the dark matter 
(r TS ). Middle and lower panels show the ratio of halo to dark matter cross-correlation coefficients without and with shot noise 
subtraction, respectively. 



of the C ; hh spectra, leads to a significant discrepancy 
between r Th and r TS . The discrepancy is w 10% at 
I = 10 for the lowest redshift cluster-sized halo sample 
(Bin 2). For the group-scale haloes (Bin 1), the deviation 
is smaller, being « 10% at I = 50, for the same redshif 
shell. However, if one subtracts shot noise from the halo 
auto-spectra (bottom panel of the figures), P s hot = 1/^h, 
then these effects can be mitigated, and the ratio r Th / r TS 
is brought within 5 < % of unity. A note of caution, is 
that we found that using the standard P s hot = to 
correct for the shot-noise lead to negative power spectra 
at high k. Since this is forbidden, we believe that such 
simple corrections are in fact an over- correction and new 
more accurate methods for accounting for the discrete- 
ness will be required [for a deeper discussion of this issue 
see [53. 



We thus conclude that the relation r Th « r TS holds to 
within 5% for I < 50, for the halo samples considerd in 
this study. This comes under the provision that the shot 
noise is accounted for and the shells are narrow. 



IX. CONCLUSIONS 



In this paper we have investigated the impact of the 
nonlinear evolution of the time rate of change of the grav- 
itational potentials on the CMB temperature auto-power 
spectrum, and also on the cross-correlation of biased den- 
sity tracers and the CMB. Linear perturbation theory 
informs us that, for nearly the entire history of the Uni- 
verse, gravitational potentials are constant and there is 
no net heating or cooling of the primordial CMB photons. 
However, at late times in the LCDM model the symme- 
try between the growth rate of density perturbations and 
the expansion rate is broken. The growth slows, and po- 
tentials begin to decay. Using the zHORIZON simulations, 
a large ensemble of iV-body simulations and analytic per- 
turbation theory methods, we explored how this picture 
changed. 

In i jllll we generated maps of the rate of change of the 
gravitational potentials at different stages in the simula- 
tion. We showed that, at redshifts z ~ 15 — 10, whilst 
the ISW signal is vanishingly small, the potentials are 
indeed evolving nonlinearly on small scales giving rise to 
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the Rees-Sciama and Birkinshaw-Gull effects - nonlinear 
infall and mass motion across the line of sight. How- 
ever, the amplitude of these effects, at these redshifts, is 
too low to be detected directly in the CMB or through a 
cross-correlation analysis. We then showed that at later 
times z > 3 the potential evolution becomes dominated 
by the large-scale ISW effect. 

In i HVI we focused on investigating the impact on the 
CMB temperature power spectrum. The late time ISW 
effect can be quantified through a line-of-sight integral 
over three power spectra: the auto-spectra of density and 
momentum, and their cross-spectrum. We used the non- 
linear PT to derive explicit expressions for each of these 
quantities. Estimates were then measured from the en- 
semble of simulations over the range z = 1-0. In all 
cases there was evidence for large-scale nonlinearity, the 
effects being strongest for the momentum auto-spectra 
and at the lowest redshifts. However, when the spectra 
were combined to produce the $ spectrum, the nonlinear 
corrections to linear predictions increased with increasing 
redshift. This was attributed to the fact that the ISW 
vanishes at early times, so leaving only the RS and BG 
effects. The standard PT was able to reproduce the non- 
linear behavior at high precision over this redshift range. 

In JV]we estimated the CMB spectrum using the Lim- 
ber approximation, we found that the nonlinear amplifi- 
cation of the ISW effect was < 10% of the linear theory 
on scales I < 50, and was also swamped by the cosmic 
variance of the linear ISW effect on these scales. On 
smaller scales the effect was more significant, however 
the primary CMB signal is more than ~ 10 3 times larger 
at this scale. We conclude that for the standard LCDM 
model, it is highly unlikely that the nonlinear ISW effects 
could contaminate the / < 1500 multipoles of the CMB 
spectrum in any traceable way. Our results support con- 
clusions from earlier studies [Tl], Qj], [zl] • 

In i jVIl we analyzed the cross-correlation of ISW with 
the dark matter density field, showing that while the 
nonlinear effects suppress this cross-correlation at early 
times, they may enhance it at very late times. This is fur- 
ther investigated in S jVIIi where we computed the ISW 
signal obtained from the cross-correlation of the CMB 
with a set of biased tracers of the density field. We 
modeled the bias using a time dependent linear model 
and also a time- and scale-dependent nonlinear model 
[H, |H, [7(| • For the biased samples we took the haloes 
measured in the simulations between z — and 1, with 
masses M > 1O 13 /i -1 M . These were then sub-divided 
into a high- and low-mass sample. The linear and nonlin- 
ear bias parameters were then estimated from the halo- 
mass cross-power spectra. The angular power spectrum 
of the ISW depends on two spectra: the cross-power spec- 
trum of the biased tracer with the mass density and the 
momentum. These spectra were estimated from the sim- 
ulations. Again there was evidence for large-scale non- 
linearity, the effects being strongest for the momentum 
cross-spectrum and at late times. The predictions from 
the nonlinear analytic PT model were found to qualita- 



tively reproduce the power spectra. On combining the 
two spectra to produce the ISW-density tracer cross- 
spectrum, we again found evidence of nonlinearity, and 
as for the case of the ISW auto-spectrum, the effects were 
more noticeable at higher redshifts. We also found that at 
late times there was an amplification of the cross-power 
spectrum. Thus at late times in the LCDM model, non- 
linear evolution can lead to a small increase in the decay 
rate of the gravitational potentials. 

In Willi we computed the angular power spectra, aver- 
aging over the halo spectra at various redshifts. We found 
that on scales I < 100 the departures from linear theory 
predictions were < 10%, and these were characterized by 
a small amplification of the signal, followed by a strong 
suppression. The departures are sub-dominant to the 
cosmic variance. We then investigated the S/J\f for the 
haloes and found good agreement with the linear theory 
expectation: the presence of bias effectively cancels out 
in the S/M expression and leads to negligible changes in 
the cross-correlation detectability. We also showed that 
through the increased Poisson noise of the biased sam- 
ple, there was a reduction in the S/Af, relative to that 
for the mass. Our analyses also demonstrated that the 
S/ N of the ISW-large-scale structure cross correlation 
is localized to a narrow angular range: more than ~ 90% 
of the overall significance arises from I < 50, or angu- 
lar scales larger than ~ 4 degrees. We therefore conclude 
that the current power spectrum analyses of Ho et al. [27| 
and Giannantonio et al. [2t| are not affected by nonlin- 
ear density evolution or scale-dependent bias to influence 
the detectability of the ISW-LSS cross-correlation. Since 
we do not repeat the exact analysis of Granett et al. [28[ 
we cannot directly address whether that result can be 
explained by nonlinear effects or whether it requires an 
alternative explanation. 

Finally, we compared the cross-correlation coefficient 
of the biased density tracers and the CMB with that 
of the dark matter and the CMB. We found that the 
relation r Th w r TS holds to within 5% for Z < 50, for the 
halo samples considerd in this study. This comes under 
the provision that the shot noise is accounted for and the 
shells are narrow. Otherwise the deviations can be large. 

The power spectrum anslysis of ISW, therefore, ap- 
pears to be a probe relatively free from contamination 
by the pernicious effects of late-time nonlinear evolution 
of the large-scale structures or scale dependent bias, at 
least for I < 100 where most of the signal is. It therefore 
continues to be a useful probe for the presence of Dark 
Energy or its alternatives [79| . 
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APPENDIX A: VALIDITY OF THE LIMBER 
APPROXIMATION 



The Limber approximation is motivated by: 

r dk e.niknlnikrz) = ^ (ri ~ r2) , (Al) 

Jo z r i 
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where the symbol S D denotes the Dirac delta function. 
Under the assumption that the spherical Bessel functions 
ji(kr) are rapidly oscillating for high enough l-s, then one 
can write an integral over a generic power spectrum as 



J 7 {l,a,ri,r 2 ) = J dk k 2 ji(kri)ji(kr 2 ) P(k) 



where the power spectrum P{k) is assumed, in a cosmo- 
logical context, to be a power law times some transfer 
function |T(fc)| 2 , P(k) = k a \T{k)\ 2 . If seen as a four di- 
mensional matrix with indices running on {7, a, ri, r2}, 
the deviation of T from a diagonal matrix in the last 
two indices, may be viewed as a measure of the error 
introduced by the Limber approximation. 

In Fig. [TS] we examine T for the case where we have 
fixed T\ to be the comoving distance to z — 0.4 and 
where r 2 varies on the X-axis. We consider three cases 
for the multipole number: I — {4,38,103}; and three 
cases for the spectral index: a — {1,-1, —3} which 
may be thought of as Pss, P^g, and P$$. We take 
fcmin = 10 _5 /iMpc _1 and fc max = 1/iMpc -1 . For the 
sake of clarity, the elements of T have been normalized 
by the maximum value of each row. Black points denote 
positive values and green ones negative entries. The di- 
agonal term (at zq = 0.4) has been marked by a vertical 
dashed line. From the figure it is clear that the devi- 
ation from a diagonal matrix is more apparent at low 
multipoles, and for more negative values of a. At higher 
I, however, the width of the T matrix shrinks around 
Zq, making the Limber approximation more precise. The 
actual error on these multipoles is related to how the 
off-diagonal terms are weighted by the time dependent 
factors, and how their sum cancels within the integra- 
tion range. 

Fig. [12] presents the errors on Cf s and Cf s , at three 
different redshifts zo — {0.2,0.4,0.6} and for thin (green 
color) and thick (red color) redshift shells (these are dis- 
played in the bottom panels). In all cases, for I > 20, the 
errors are below 3%. We find that for Cf s , the net result- 
ing error is larger for thin redshift shells than thick ones. 
This is inverted for Cf 5 , where the contribution to the 
off-diagonal terms are smaller for thin shells. However 
errors remain always below the few-percent level. The 
amplitude of the errors are defined by: the actual width 
of the peak around z = zq; the amplitude of the oscillat- 
ing floor around the wings of the peak at z = z$\ and the 
actual width of the redshift integration range compared 
to the width of the peak at z = zq. 
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FIG. 18: Rows of the matrix IF (I, a, n, ra) (normalized by the diagonal term) corresponding to n (zo = 0.4) versus the redshift 
corresponding to ra, under different choices of I and a. Given the logarithmic scale, green color displays negative values, black 
points positive ones. 
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FIG. 19: Comparison of the exact Ci evaluation for the ISW-density tracer correlation with the Limber approximation evalu- 
ation. The top three panels show the relative errors for a near, intermediate and far density tracer survey. Solid lines denote 
predictions for thick redshift shell and dash lines denote thin redshift shells. The corresponding bottom panels show the redshift 
distributions. 



